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Abstract 

We present an application of equation-free computation to the coarse-grained feedback linearization 
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illustrative example is a kinetic Monte Carlo realization of a simplified heterogeneous catalytic reaction 
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1 Introduction 



A fundamental prerequisite for the design of control systems is the availability of reasonably 
accurate closed form dynamical models. Typically, such models arise in the form of evo- 
lution equations (ordinary differential, differential algebraic, partial differential, possibly 
integrodifferential equations). Such equations are typically derived from conservation laws 
(e.g. mass, momentum and energy balances) closed through constitutive equations (e.g. 
Newtonian stresses in fluid flow, or mass-action kinetics expressions for chemical reactions) ; 
system identiflcation may also play a role in obtaining and/or closing such macroscopic 
models. Many real-world problems of current engineering interest are characterized -due 
to their stochastic/microscopic nature, and nonlinear complexity- by the lack of such good 
explicit, coarse-grained macroscopic evolution equations. Instead, the underlying physics 
description may be available at a much finer, more detailed level: the evolution rules may 
be given in the form of molecular dynamics, kinetic Monte Carlo, Markov-chain or hybrid 
schemes. When this is the case, conventional continuum algorithms cannot be used directly 
for systems level analysis and controller design. Bridging systematically the enormous gap 
between microscopic space and time scales of a complex physical/material system descrip- 
tion and the macroscopic ones at which we want to design and control its behavior is a 
grand challenge for modeling and computation. Over the past few years we have demon- 
strated that an equation-free approach (based on coarse timesteppers) [Thcodoropoulos et 
al., 2002; Makeev et al.,2002; Kevrekidis et al., 2003; Siettos et al., 2003b; Kevrekidis et 
al., 2004], can establish a link between traditional continuum numerical analysis and micro- 
scopic/ stochastic simulation. This is a mathematics-assisted computational methodology, 
inspired from continuum numerical analysis, system identiflcation and large scale itera- 
tive linear algebra, which enables microscopic-level codes to perform system-level analysis 
directly, without the need to pass through an intermediate, coarse-grained, macroscopic- 
level, "conventional" description of the system dynamics. The backbone of the method 
is the on-demand identiflcation of the quantities required for continuum numerics (coarse 
residuals, the action of coarse slow Jacobians, eigenvalues, Hessians, etc). These are ob- 
tained by repeated, appropriately initialized calls to an existing fine scale time-stepping 
routine, which is treated as an input-output black box. The key assumption is that de- 
terministic, macroscopic, coarse models exist and close for the expected behavior of a few 
macroscopic system observables, yet they are unavailable in closed form. These observables 
(coarse-grained variables) are typically a few low moments of microscopically evolving dis- 
tributions (e.g. surface coverages, the zeroth moments of species distributions on a lattice 
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model of a surface reaction). 

The present work aims at developing a systematic approach to the feedback regulator 
synthesis problem, where both the closed-loop dynamics linearization and pole-placement 
objectives are simultaneously attained by using the equation-free timestepper methodology. 
The feedback linearization and the pole-placement objectives for the unavailable coarse- 
grained dynamics are met in a single-step, circumventing the lack of an explicit dynamic 
process model. The proposed approach is illustrated through the use of a coarse time- 
stepper based on a kinetic Monte Carlo realization of a simplified surface reaction scheme 
for the dynamics of NO oxidation by H2 on Pt and Rh surfaces. The present paper 
is organized as follows: In section 2 we briefly discuss the traditional nonlinear control 
methodologies that rely on the notion of feedback linearization along with the associated 
restrictions encountered at the implementation stage. In section 3 we succinctly review a 
recently proposed approach that allows the attainment of both the feedback linearization 
and pole placement objectives in a single step, effectively overcoming the restrictive condi- 
tions associated with the classical exact feedback linearization approach. In section 4 the 
interplay of the proposed nonlinear control procedure with coarse timesteppers is outlined, 
and the natural integration of the respective frameworks illustrated. Section 5 presents the 
simulation results using the proposed methodology on an illustrative kinetic Monte Carlo 
model, followed by some concluding remarks in section 6. 

2 Fundamentals of feedback linearization of nonlinear systems 

In order to meet a set of performance specifications or design objectives, process control 
introduces feedback to appropriately modify the dynamics of a system. Placing the closed- 
loop poles at desirable locations in the complex plane, and thus shaping the closed loop 
system dynamics and time constants, is a popular controller synthesis method for linear 
systems, in part, due to its intuitive appeal [Chen, 1984]. Typically one requires fast decay 
of the closed loop variables to their nominal steady state values; yet the design should 
not lead to high feedback gains due to possible saturation problems. Fine-tuning of the 
closed-loop eigenvalues is performed in practice through a combination of optimization 
techniques, heuristic rules and trial-and-error approaches [Chen, 1984]. Traditional pole- 
placement state feedback control laws for nonlinear systems are based on local linearization 
around a reference steady state, and the subsequent use of linear design methods. The re- 
sults are, of course, only locally valid, and may lead to unacceptable performance, even in 
the presence of only mild nonlinearities. Nonlinear feedback control laws thus need to be 
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derived, capable of directly coping with the system nonlinearities. A pole-placing feedback 
regulator should be capable of bringing the system/process state back to the design steady 
state in a fast and smooth manner in the presence of disturbances; if the design steady 
state is unstable, the primary control objective is its stabilization. In the pertinent body 
of literature two main model-based pole-placing controller synthesis methods emerge, both 
based on geometric control theory. The first one is exact input/output (I/O) feedback lin- 
earization, where the introduction of nonlinear state feedback induces linear I/O behavior 
of the system of interest, forcing the system's output to follow a prespecified linear and 
stable trajectory. This approach directly generalizes the linear result of placing the closed- 
loop poles at the system's zeros and at a set of prespecified values, and is restricted within 
the class of minimum-phase systems [Isidori, 1999]. Regulation and/or stabilization of a 
system/process, however, is understood in terms of forcing the system's state to return to 
the design steady state (if driven away from it in the presence of disturbances). Further- 
more, process output tracking problems for step changes in the output set-point values, 
can be easily reformulated as regulation problems relative to the equilibrium point that 
corresponds to the final set-point value. The second approach is geometric exact feedback 
linearization, traditionally implemented in a two-step design procedure [Isidori, 1999]: A 
simultaneous implementation of a nonlinear coordinate transformation and a state feedback 
control law in the first step transforms the original nonlinear system to a linear and con- 
trollable one. Well-established finear pole-placement techniques for the transformed finear 
system can be used in the second step. However, the aforementioned classical geometric 
exact feedback linearization approach relies on a set of very restrictive conditions, that can 
hardly be met by any physical system. 

In this work a systematic approach to feedback regulator synthesis is proposed for the 
coarse-grained dynamic behavior of systems described by atomistic/stochastic ("fine scale") 
simulators. The closed-loop dynamics linearization and the pole-placement objectives are 
simultaneously attained using the equation-free timestepper-based methodology. Note that 
our primary control objective is to assign the closed-loop eigenvalues rather than shaping 
the entire I/O behavior of the system under consideration. Furthermore, applying the 
methodology introduced in [Kazantzis, 2001], we investigate the possibility of circumvent- 
ing the set of restrictive conditions associated with the two-step classical exact feedback 
linearization approach, by meeting the feedback linearization and the pole-placement objec- 
tives in a single-step, and without being limited by the availability of an explicit dynamic 
process model. 
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3 Mathematical Preliminaries - Problem Formulation 



In the context of the present study, the system dynamics are described by a nonhnear 
discrete-time macroscopic ( "coarse" ) model of the form: 

x{k + l) = ^{x{k),u{k)). (1) 

Here k e — {0, 1, ...} is the discrete-time index, x{k) e i?" is the vector of (coarse) 
state variables, u{k) e i? is the manipulated input variable and ^{x, u) represents a vector 
function defined on i?" x R. In our case this function is not known, and will be identified 
on the fly with the aid of the fine scale simulator. Without loss of generality, it is assumed 
that the origin = is an equilibrium point (coarse steady state) of (1), that corresponds 
to vP — 0: $(0, 0) =0. If a non-zero coarse steady-state (x°, m°) ^ (0, 0) is located, then a 
simple transformation: x — x — x^ ,u — u — u^ will map it to the origin in the new coordinate 
system. Let F be the Jacobian matrix of m) evaluated at x = 0: F — —(0,0), and 

G the vector: G — -;^(0, 0) which is assumed to be non-zero. The following assumption is 
ou 

also made: 

Assumption I: The nx n matrix: 

C=[G\FG\...\F''-^G] (2) 

has rank n. This implies that the coarse linearization of (1) around the origin a; = is 
controllable [Isidori, 1999]. 

It is appropriate, at this point, to briefiy review and outline basic features of the classical 
exact feedback linearization approach in the discrete-time domain. In the first step, and 
under a set of rather restrictive conditions [Aranda-Bricaire et al., 1996; Califano et al., 
1999; Grizzle, 1986; Jacubczyck, 1987; Lee et al, 1987; Lin and Brynes, 1995; Nam, 1989], 
a nonlinear coordinate transformation: z — T{x) is sought along with a state feedback 
control law: u — ^{x,v) (with v being an external reference input), such that the original 
system (1) is transformed to the following linear one: 

z{k + 1) = Az{k) + hv{k) (3) 

where {A, h) is a Brunowsky controllable pair of matrices [Chen, 1984]. In the second step, 
standard linear pole-placing feedback techniques are used to arbitrarily assign the poles 
(equivalently the time-constants) of the closed-loop system. In particular, a constant-gain 
vector K is calculated, such that the state feedback law: v — —Kz induces the desirable 
closed-loop dynamics: 

z{k + l)^Az{k) = {A-hK)z{k) (4) 
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with A = A — bK being the closed-loop system's characteristic matrix with prescribed 
eigenvalues. 

At this point it would be appropriate to review an alternative single-step design method 
for linear systems: 

x{k + l) = Ax{k) + hu{k), (5) 

where A, b are constant matrices with appropriate dimensions, that was first introduced by 
D. Luenberger (1963). This alternative approach serves as the methodological basis for the 
development of a nonlinear analogue introduced in [Kazantzis, 2001] and briefly outlined 
in the next section. According to the ideas reported in [Luenberger, 1963] a single-step 
simultaneous implementation of a linear coordinate transformation: z — Tx coupled with 
a hnear state feedback control law: u — —Kz is sought, that induce the following closed- 
loop dynamics: 

z{k + 1) = Az{k) (6) 

A is the closed-loop system's characteristic matrix that carries the prescribed set of eigen- 
values due to the control law applied. This requirement can be mathematically translated 
into a quadratic matrix equation that the unknown transformation matrix T should satisfy: 

TA-AT^ TbKT (7) 

If T is non-singular (invertible), one can easily show that the inverse transformation matrix 
W — satisfies the following linear matrix equation: 

AW-WA^ bK. (8) 

It is known from linear algebra that, if matrices A and A have disjoint eigenspectra, the 
above matrix equation (8) admits a unique solution W [Chen, 1984; Gantmacher, I960]. 
Furthermore, invertibility of the solution can be ensured iff the pair of matrices (^4, b) is 
controllable and the pair {K, A) is an observable one [Chen, 1984]. As shown in [Luenberger, 
1963], if T is the unique invertible solution to the matrix equation (7), then the linear state 
feedback control law expressed in the original variables x 

u{k) = -KTx{k) (9) 

induces the closed-loop dynamics: 

x{k + 1) = T-^ATx{k) (10) 

Since matrices T~^AT and A are similar, it can be easily inferred that the closed-loop 
system has the desirable set of poles assigned by the control law (9). 
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3.1 Single-Step Feedback Linectrization With Pole- Placement 

Motivated by D. Luenberger's linear approach [Luenberger, 1963], let us now succinctly 
review the ideas presented in [Kazantzis, 2001] on its nonlinear generalization. One seeks 
to simultaneously implement a nonlinear coordinate transformation, z — S{x) coupled with 
a nonlinear state feedback control law, u — —cz — —cS{x), where c is an arbitrary constant 
row vector (a design parameter of the proposed method) that induce linear closed-loop z- 
dynamics: 

z{k + l)^Az{k). (11) 

The poles of the closed-loop dynamics (11) are realized by the eigenvalues of the arbitrar- 
ily prescribed matrix A: the characteristic matrix of the linear closed-loop dynamics (11). 
Therefore, the eigenspectrum of A should be judiciously selected to favorably shape the 
dynamic characteristics of the controlled system's response. In the nonlinear case, these de- 
sign requirements arc embodied into the following system of nonlinear functional equations 
(NFEs) that need to be satisfied by the unknown transformation map S{x): 

Smx,-cS{x))) = AS{x) 

S{0) = 0. (12) 

The accompanying initial condition S{0) — merely reflects the fact that equilibrium 
properties must be preserved under the proposed coordinate transformation. 

For the study of the mathematical properties of the solution of the NFEs (12) and 
within the class of real analytic systems, a number of assumptions are essential as shown 
in [Kazantzis, 2001]: 

Assumption II: All the eigenvalues ki,{i — 1, ...,n) of matrix A should he inside the 
unit disc on the complex plane (stability requirement imposed on the closed-loop dynamics 

(11))- 

Assumption III: The eigenspectra a{A), a{F) of matrices A and F respectively should 

be disjoint: a{A) n a{F) = 0. 

Assumption IV: The eigenvalues ki of A should not be related to the eigenvalues 
Xj,{j = 1, n) of F through any equations of the type: 

n = A, (13) 

i=l 

{j — 1, ...,n), where all the m^'s are non- negative integers that satisfy the condition: 

n 

1=1 
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Assumption V: The pair of matrices (c, A) is chosen such that the following matrix 

O: 

(15) 

has rank n: rank(O) = n (Observability condition imposed on (c, A)). 

Lemma: [Kazantzis, 2001] For a real analytic system (1), let the above assumptions I-V 
hold true. Then, the system of NFEs (12) admits a unique locally analytic and invertible 
solution z — S{x) in a neighborhood of the origin x — 0. 

We include here a number of remarks discussing the conditions and imphcations of this 
Lemma; a more detailed discussion can be found in [Kazantzis, 2001]. 

Remark 1: The "non-resonance" conditions (13) and (14) are required for the existence 
of a unique formal powcr-scrics solution to the system of NFEs (12). The assumption for 
the eigenspectrum of matrix A to lie inside the unit disc plays a key role in the uniform 
convergence of this formal power-series solution in the neighborhood of the origin x = 
with a non-zero radius of convergence, and thus for the solution's analyticity. Finally, 
Assumptions I and V are necessary and sufficient conditions for local invertibility of the 
solution. 

Remark 2: It is useful to consider the linear case: $(x, u) = Fx + Gu where F, G are a 
constant matrix and vector of appropriate dimensions respectively. In this case, the unique 
solution of the system of NFEs (12) is, w — Sx, where S is the solution to the quadratic 
matrix equation: 

SF-AS^ SGcS. (16) 

which coincides with (7) in D. Luenberger's analysis. Please notice, that under the as- 
sumptions stated the above matrix equation (16) admits a unique and invertible solution 
S [Chen, 1984]. 

Let us now consider: z — S{x) to be the solution to the associated system of NFE's 
(12) defined in a neighborhood of x = 0. It has been shown in Kazantzis (2001) that the 
simultaneous implementation of the nonlinear coordinate transformation: z = S{x) and 
the nonlinear state feedback control law: 

u{k) = -cS{x{k)) (17) 

results in linear closed-loop ^-dynamics: 

z(k + 1) = Az{k) (18) 



O 



c 

cA 



cA' 
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whose poles are realized by the eigenvalues of matrix A. Indeed, one can easily show that 
the closed-loop system dynamics expressed in the 2;-coordinates satisfy: 

z{k + l) = S{x{k + 1))^ Smx{k),-cS{x{k)) 

= AS{x{k)) = Az{k). (19) 

Remcirk 3: Note that in the linear case, one calculates a feedback control law: 

u{k) = -cSx{k) (20) 

where S is the solution to (16), that induces the following closed-loop dynamics: 

x{k + 1) = (F - GcS)x{k) (21) 

Using equation (16), the closed-loop dynamics (21) can be rewritten as follows: 

x{k + 1) = {S-^AS)x{k) = Ax{k). (22) 

Notice that A,A = S'^AS are similar matrices, and therefore, the closed-loop system (22) 
has the desirable poles. One can consider this approach as the natural extension of D. 
Lucnbcrger's linear result for pole-placement (10) to nonlinear systems. 

Remark 4: The graph of the mapping z = S{x) is rendered invariant for the com- 
posite system (1) and (11) under the state feedback control law: u{k) = —cS{x{k)) [Carr, 
1981]. Furthermore, the system of NFEs (12) represent the associated invariance functional 
equations for the composite system (l)-(ll) [Guckenheimcr and Holmes, 1983], and the re- 
striction of the composite system dynamics under the above feedback law on the invariant 
manifold/ solution of (12) coincides with the linear closed-loop dynamics (11). 

Remark 5: The primary idea of the proposed single-step design approach is to avoid 
the intermediate step of transforming the original system into a linear controllable one 
with an external reference input, which allowed us to circumvent the restrictive conditions 
associated with the classical exact feedback linearization method [Aranda-Bricaire, 1996; 
Califano, 1999; Grizzle, 1986; Jakubczyck, 1987; Lee, 1987; Lin and Byrnes, 1995; Nam, 
1989]. It should be pointed out, that the design method described does not involve an 
external reference input, however, and therefore other control objectives such as trajectory 
tracking or model matching can not be met [Isidori, 1999]. 

In the present study the NFEs (12) will be solved using the equation-free computational 
framework. However, for completeness and comparative accuracy, one needs to employ a 
alternative solution scheme/method for the system of NFE's (12). This method involves 
expanding $(x, u) as well as the unknown solution S{x) in a Taylor series and equating 
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the Taylor coefficients of tlie same order of botli sides of the NFE's (12). This procedure 
leads to linear recursion formulas, through which one can calculate the A^-th order Taylor 
coefficients of S{x), given the Taylor coefficients of S{x) up to the order — 1. As shown 
in [Kazantzis, 2001], in the derivation of the recursion formulas, it is convenient to use the 
following tensorial notation: 

a) The entries of a constant matrix A are represented as aj , where the subscript i refers 
to the corresponding row and the superscript j to the corresponding column of the matrix. 

b) The partial derivatives of the /U-th component ^^{x,u) of the vector function ^{x,u) 
with respect to the state variables x evaluated at {x, u) — (0, 0) are denoted as follows: 

etc., where z, j, k, ..=1, n 

c) The partial derivatives of the ^-th component ^^{x,u) of the vector function (^[x,u) 
with respect to the input variable u evaluated at (x, u) — (0, 0) are denoted as follows: 

(24) 

etc. 

d) The standard summation convention where repeated upper and lower tensorial indices 
are summed up. 

Under the above notation the l-th. component Si{x) of the unknown solution S{x) can 
be expanded in a multivariate Taylor series as follows: 

1 

N\ 



+ ^ Xi^Xi^...Xij^ + ... (25) 



Similarly one expands the components of the vector function $(x, u) in multivariate Taylor 
series. Substituting the Taylor expansions of S{x) and ^{x,u) into (12) and matching the 
Taylor coefficients of the same order, the following relation for the AT-th order terms may 
be obtained [Kazantzis, 2001]: 

N 

E E ^r-^M/r-/r-^r-</)=<^^-^" (26) 

L=\ 0<.m-i <'m2^---^iTij^ 
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where: 

<^ = E E g^ysr-'''' (27) 

P=l 0<n]^ <'>T'2<-<np 

ii,...,ijs[ = l,...,n and / = l,...,n. Notice that the second summation symbol in (26) 

AT! 

(and similarly in (27)) suggests summing up the relevant quantities over the — ; r 

mi\...mLi 

possible combinations to assign the A'^ indices (ii, ijv) as upper indices to the L positions 
{/jj, ■■■fjj^} (and {tTjj, •■•tTj^}), with mi of them being put in the first position, m2 of them in 

L 

the second position , etc. (^rrii — N). Moreover, notice that equations (26,27) represent 

i=l 

a set of linear algebraic equations in the unknown coefficients S^^''"'^^ for N > 2. For 
N = 1, equations (26,27) yield the quadratic matrix equation (16) (or (7) in D. Luenberger's 
approach). It should be pointed out, that the above scries solution method for the NFE's 
(12) is amenable to a computer-based implementation and can be readily carried out in an 
automatic fashion with the aid of a symbolic software package such as MAPLE. 

4 An Equation- Free Approach to the Feedback Linearization Prob- 
lem 

As shown in the previous section, the system of NFEs (12) admits a unique analytic solu- 
tion. However, such an analytic transformation is difficult to derive in the general case, and 
a numerical solution scheme becomes necessary. We will now assume that the model equa- 
tions are not explicitly available, but we do have a "black box" subroutine that, given the 
state of the system G i?" , "Uq € -R at time tk = kT reports the result of the system after a 
time horizon T (i.e., will report x{tk+i = {k + l)T) = ^t{xo,Uo)). This subroutine could be 
a "legacy" dynamic simulator; alternatively, it can be a "coarse timestcppcr" involving the 
lift, run and restrict steps discussed briefly below and in more detail in [Makeev ct al., 2002; 
Gear et al., 2002; Kevrekidis et al. 2003, Siettos et al., 2003b]. The coarse timestepper, 
which we use in the equation-free framework for coarse-grained controller design (for linear 
quadratic control, pole placement and feedback linearization) [Siettos et al., 2003a, Siettos, 
et al., 2004a, Armaou et al., 2004a, 2004b, Siettos et al., 2004b] consists of the following 
elements (Figure 1): 

• a lifting operator /i, transforming a macroscopic initial condition (typically zeroth- 
or first-order moments of the microscopically evolving distributions) to one (or more) 
consistent microscopic realizations; 

• evolution of the microscopic realizations using the microscopic simulator for an appro- 
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priately chosen (relatively short) macroscopic time T, {the reporting horizon). 

• a restriction operator M, transforming the resulting microscopic distributions back to 
the macroscopic description (obtaining their macroscopic observables) . Lifting from 
microscopic to the macroscopic and then restricting again should have theoretically 
no effect (modulo roundoff), that is, fj,M — I. 

This coarse timcstcppcr, appropriately initialized and executed can serve in the "on 
demand" estimation of model right-hand-sides, of the action of "coarse slow" Jacobians as 
well as derivatives with respect to parameters, in the computation of coarse fixed points 
and their leading eigenvalues - in short, of exactly the quantities that a linear or nonlinear 
controller design algorithm would need evaluated through a macroscopic model (had such 
a model been available) to perform its task. 

For our problem, we use the coarse timestepper in a coarse fixed point algorithm to 
converge on the desired coarse nominal equilibrium xq; we then proceed as follows (re- 
markably, the algorithm is the same for the case of legacy dynamic simulators and coarse 
timesteppers of microscopic/stochastic models): 

• Discrctizc the domain C R"- of the state-space, where a numerical solution of the 
feedback linearization problem is sought in a mesh of, say, N points. 

• Write the transformation vector S{x) as a power series expansion up to order p around 
the equilibrium xo i.e. write S{x) as S{x; h), where h e R"^ is the vector of the power 
series coefficients. For example for a 2-dimensional problem S{x; h) = S{xi, X2; h) can 
be written as: 

Si{xi,X2; ai=i,...p) = aixi + a2X2 + -^asxl + ^04X3 + 05X1X2 + ... + 0{p + 1) (28) 
S'2(xi, X2; 6j=i,...p) = 61X1 + 62X2 + -^hxl + ^bixl + 65X1X2 + ... + 0{p + 1) (29) 
where: 

h = [oi, 02, ttp, bi, &2, •-, bp] (30) 
Then, write the feedback control law as in (17). 

• Calculate the values of the unknown coefficients of ^(x; h) using a matrix-free iterative 
nonlinear solver [Kelley, 1999], or possibly an unconstrained optimization algorithm, 
such as the Broyden, Fletcher, Goldfarb, Shanno (BFGS) method. 
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The optimization problem can be stated as finding the values of the vector h such that 
the sum of squared errors on the discretization mesh is minimized, i.e.: 

1 ^ 

min-Y. \\G,{h) \\l (31) 

where the vector function Gi{h) is defined as: Gi{h) — S{^T{xii —cS{xi); h) — AS{xi, h), 
\/xi on the discretized mesh, and || • II2 is the standard Euclidean norm in the above 
minimum norm problem [Luenberger, 1969]. 

The quantities involved in the optimization computations (e.g. the values Gi) are eval- 
uated repeatedly using the (legacy or coarse) timestepper for each value Xi in the mesh. 
Remcirk 6: The single-step feedback linearization problem under consideration admits an 
alternative formulation, where the inverse transformation map: x — w{z) is sought that 
satisfies the following system of NFEs: 

w{Az) = ^{w{z),-cz) 
w(0) = (32) 

where: 

X = w{z) = S-\z) (33) 

The above functional equation is structurally simpler (first-order) than (12) (second-order), 
since in the latter the unknown vector function S{x) appears through two consecutive 
function composition operations. Furthermore, it can be easily shown that the above 
problem reformulation leads to the same results, namely the same feedback linearizing 
control law [Kazantzis, 2001]. Notice, that in this case we expand w{z) (instead of S{x)) 
in a power series and we then seek the values of the vector h' such that the sum of squared 
errors on the discretized domain (w.r.t z state-space, say D'" C R^) is minimized, i.e.: 

1 ^ 

min-El|G':(/i')||^ (34) 
^ i=i 

where the vector function G[{h') is defined as: G[{h') — w{Azi) — ^T{w{Azi), —czi); h'), 
\/zi on the discretized mesh. 

Upon convergence we find the desired transformation S{x) symbolically by applying a 
functional inverse on w{z). More generally, matrix-free iterative linear algebra approaches 
can be used to solve the discretized nonlinear functional equations; in these methods the 
action of the Jacobian is estimated by appropriately initialized nearby initial conditions 
(dictated, for example, by a GMRES protocol). It is worth noting that, if the problem 
dynamics are characterized by a separation of time scales, and the long-term dynamics lie on 
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a low dimensional, attracting manifold, the dynamical integration involved in timestepping 
may be beneficial to the convergence of such iterative solution techniques [Wington et al., 
1985; Kelley et al., 2004]. 

5 An Illustrative Case Study 

5.1 The Deterministic Version 

Our illustrative example consists of a simplified mechanism for the dynamics of NO reduc- 
tion by H2 on Pt and Rh surfaces. The simplified deterministic mean field model for this 
mechanism is given by: 

— = a{\ — x) — 'fx — u{l — x^x = L{x, u) (35) 

where x is the coverage of adsorbed NO, a is the rate constant for NO adsorption (incor- 
porating the gas phase NO partial pressure), 7 is the rate constant for NO desorption, 
and u is the reaction rate constant (and, in our scheme, the control variable). In order to 
transform the problem back to the discrete time formulation, we take a forward Euler step 
of the continuous time problem 

x{k + 1) = x{k) + TL{x{k),u{k)) = ^{x{k),u{k)). (36) 

Simulation results were obtained for: a — l,j — 0.01. This model, exhibits two regular 
turning points (at u ~ 3.96 and u ~ 26) as shown in the bifurcation diagram (Figure 2). 
We want to derive a nonlinear feedback control based on the proposed methodology, to 
stabilize the timestepper at the open- loop unstable stationary state {xq = 0.5559, = 4). 

We chose T = 0.1 as the reporting time horizon; the open loop eigenvalue at the nominal 
steady state is 0.1459, and the characteristic time is 6.85; A time step of 0.1 is therefore 
sufficient for accuracy of the Euler integration step and numerical stability. 

5.2 The Microscopic/Stochastic Version 

The procedure remains essentially the same when the timestepper results are obtained 
through short bursts of microscopic simulation. Here for the stochastic simulations of the 
mechanism embodied in (36) we used the Gillespie Stochastic Simulation Algorithm (SSA) 
[Gillespie, 1976, 1977]. 

Given the value of the surface coverage at time t — we computed the expected value 
of the coverage after a reporting time horizon T by simulating a system with a relatively 
large number of available sites (say Nsize), averaging over several realizations (say Nmn)] the 
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system size and number of realizations were chosen here to be Ngi^e = 100^ and A^run = 100, 
respectively. The time horizon was again selected to be T = 0.1. The Monte Carlo model 
is considered as a "black box" coarse timestepper x{k + 1) = $7-(x(A;), ^(A;)). 

The coarse Jacobian (here, a single derivative, which doubles as the coarse eigenvalue) 
at the fixed point is estimated by wrapping a Newton's method around the coarse KMC 
timestepper. The coarse identified model (Jacobian and right hand-side) is then used for 
tracing the solution branch by coupling to a pseudo-arc-length continuation scheme [Keller, 
1977]. For the continuation we used Ngize — 200^ and A^^un = 1000. For details on the 
computation of coarse stationary states and coarse bifurcation diagrams in an equation- 
free framework see [Makeev et al., 2002; Gear et al., 2002; Kevrekidis et al. 2003]. The 
resulting bifurcation diagram coincides with the one obtained through the deterministic 
timestepper. Given the unstable coarse stationary state at w = 4, the requisite functional 
equation for simultaneous feedback linearization and pole placement was solved using the 
coarse timestepper and minimizing Equation (30) using the BFGS method. To implement 
this procedure we used deviation variables defined and u' — u — Uq, while A 

now is a scalar chosen as 0.8. 

We derived the unknown transformation map S{x) numerically by the following distinct 
ways: 

a) Analytically, by expanding S{x) in a power series and retaining quadratic terms of the 
form S{x) — aiX-\-Q.ha2x'^ , substituting u — —S{x) into ^{x{k),u{k)) and then expanding 
^{x{k),—S{x)) in Taylor series around the equihbrium (0,0). The values of the unknown 
coefficients a^'s are computed by equating terms of the same order on both sides of NFEs 
(12). 

b) Equation-free by using the "black-box" KMC timestepper approach, i.e. by solving 
the optimization problem as appearing in (31) using the BFGS quasi-Newton method and 
a line search technique. 

Here the domain of interest was chosen as D & R= [—0.1 0.1] and was discretized into 
25 equally spaced points. In Figure 3 we plot the derived S{x). 

The transformation found by solving the NFE using timestcpping is later used to close 
the loop (simultaneously linearizing and assigning poles for the closed loop system dy- 
namics) Figure 4 demonstrates responses resulting from the desired closed loop dynamics 
z{k + 1) = S{k + 1) = 0.8z{k) = 0.8S{k) (dotted lines) and that of the numerically 
obtained transformation S{x) when applying the control law on the coarse KMC timestep- 
per (solid ones). Figure 5 shows the dosed loop responses of the deterministic mean field 
model and of the Kinetic Monte Carlo version starting from different initial conditions. 
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These were obtained through the solution of the norm minimization problem (31) using 
the deterministic and the stochastic KMC model respectively. The feedback linearizing 
transformation was found by minimizing (31) using the BFGS method. The obtained re- 
sults confirm the effectiveness of the proposed equation-free nonlinear controller design 
methodology, demonstrating successful stabilization and regulation of the process at the 
unstable stationary state. 

6 Concluding Remarks 

We demonstrated how feedback linearization with pole placement in a single step, analyzed 
in [Kazantzis, 2001] for closed form equation models, can be performed in an equation-free 
framework by acting directly on a fine scale simulator. The illustrative example used a 
stochastic realization of a simplified model of a catalytic surface reaction. Admittedly, 
the example is a very simple one; in particular, it is (coarsely) one-dimensional, and for 
such systems a feedback linearization transformation always exists [Isidori, 1999]. Yet the 
timestepper based, equation-free methodology illustrated is not restricted to one dimen- 
sional (when coarse-grained) problems; all the elements of the method (the timestepper, 
the location of unstable fixed points and their leading slow eigenvalues, the solution of 
the corresponding functional equation) remain effectively the same in higher-dimensional 
cases. More sophisticated matrix-free iterative methods can be used to solve the requisite 
functional equation by acting directly on the fine-scale simulator. Parallel computation (a 
different replica fine scale simulation of the same initial condition performed on each proces- 
sor) and computational tools like In situ Adaptive Tabulation (IS AT, [Pope, 1997]) can be 
used to alleviate, when appropriate, the computational wall clock time and effort required 
to estimate the necessary coarse-grained quantities. Furthermore, if a strong separation 
of time scales (a spectral gap) appears in the coarse-grained dynamics, and the long-term 
behavior lies on a low-dimensional "slow manifold" , it is possible to take advantage of this 
through timestepping to solve an effective NFE of reduced dimension [Gear et al., 2004]. In 
this paper we assumed that we knew what "the right" macroscopic observable was, in which 
to restrict the microscopic system dynamics; the detection of appropriate such observables, 
either through data analysis [Coifman et al., 2004] or observer design is a subject we are 
currently pursuing. 
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Figure Captions 



Figure 1: Schematic of the coarse timestepper in a controller design framework 

Figure 2: (a) Coarse Bifurcation diagram of the kMC model, obtained by the coarse 
timestepper, (b) blow up of the diagram near the equilibrium of interest; solid lines corre- 
spond to stable coarse steady states while the dotted ones correspond to unstable coarse 
steady states 

Figure 3: S{x) as computed analytically (solid line) and using the black-box coarse 
KMC timestepper (dotted line) 

Figure 4: Transients of S{k + 1) = 0.8S'(A;), corresponding to the desired closed loop 
dynamics, (solid lines) and S{x{k)) using the computed control law on the coarse KMC 
timestepper (dotted lines) 

Figure 5: (a) Transient response for 0.1 initial perturbation of the coarse state variable 
from the coarse equilibrium, (b) Transient response for 0.2 initial perturbation of the coarse 
state variable from the coarse equilibrium, (c) Transient response of the control variable 
for 0.2 initial perturbation of the coarse state variable from the coarse equilibrium (lower 
ones correspond to -0.2). Simulation runs for the KMC timestepper were obtained with 
= 100^ and N„,„ = 100. 
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Figure 2: (a) Coarse Bifurcation diagram of the kMC model, obtained by the coarse timestepper, (b) blow 
up of the diagram near the equilibrium of interest; solid lines correspond to stable coarse steady states 




Figure 3: S{x) as computed analytically (solid line) and using the black-box coarse KMC timestepper 
(dotted line). 



Six) \ 




"'o 50 100 150 

time (s) X 0. 1 



Figure 4: Transients of 5'(A; + 1) = 0.8S{k), corresponding to the desired closed loop dynamics, (solid lines) 
and S{x{k)) using the computed control law on the coarse KMC timestepper (dotted lines) 




Figure 5: (a) Transient response for 0.1 initial perturbation of the coarse state variable from the coarse 
equilibrium, (b) Transient response for 0.2 initial perturbation of the coarse state variable from the 



